#include "AN_rc_w_lnu.h"
bool AN_RC_w_lnu(CDraw &para, float wmass, float mw_d, float mw_u, TLorentzVector &Vlep, MissingET* Vmet, TLorentzVector &re_w, TLorentzVector &re_nu){
	para.status.Set_Current_Position("begin AN_RC_w_lnu");
	int lnum,mnum,key;


	TLorentzVector pneutrino=TLorentzVector(Vmet->MET*std::cos(Vmet->Phi),Vmet->MET*std::sin(Vmet->Phi),0,0);
	TLorentzVector pw_choosed, pnu_choosed;



	float re_mass=0;
	//float xi_min=10000;
	bool Jre_w=false;

/// ShowMessage(2,"***************");
/// ShowMessage(2,"begin re W",wmass, mw_d, mw_u);
///	ShowMessage(2,"momentum",pneutrino, Vlep);
    int count=0;
	for(float wmass_tmp=0;wmass_tmp<5*wmass;wmass_tmp=wmass_tmp+1){
		TLorentzVector pw,  pnu;
		Jre_w=AN_RC_lnu(para, wmass_tmp, Vlep, pneutrino, pw, pnu);
////    ShowMessage(2,"wmass_tmp",wmass_tmp, Jre_w);
////    ShowMessage(2,"pnu",pnu);
////    ShowMessage(2,"pw",pw);
  		if(!Jre_w){
			continue;
  		}
		else{
	    	re_mass=pw.M();
	    	pw_choosed=pw;
	    	pnu_choosed=pnu;
			break;
		}
	////float xi=std::abs(pw.M()-wmass);
	////if(xi<xi_min){
	////	xi_min=xi;
	////	re_mass=pw.M();
	////	pw_choosed=pw;
	////	pnu_choosed=pnu;
	////}
	////else{
	////	count++;
	////}
	////ShowMessage(2,"wmass_tmp",wmass_tmp, xi, xi_min);
	////if(count>10){
	////	break;
	////}
	////std::string position="xi: "+Float_to_String(xi)+" xi_min: "+ Float_to_String(xi_min) + " re_mass: " + Float_to_String(re_mass)+ " pw.M: " + Float_to_String(pw.M());
	////para.status.Set_Current_Position(position);
  	}

	bool Jre_cut=false;

	if(re_mass>mw_d&&re_mass<mw_u){
		Jre_cut=true;
	}
	if(Jre_w&&Jre_cut){
		re_w .SetPxPyPzE(pw_choosed.Px(),pw_choosed.Py(),pw_choosed.Pz(),pw_choosed.E());
		re_nu.SetPxPyPzE(pnu_choosed.Px(),pnu_choosed.Py(),pnu_choosed.Pz(),pnu_choosed.E());
		if(re_w.M()<10){
			para.status.Set_Show_Status(true);
		}
		std::string position="has correct W  remass: "+Float_to_String(re_mass);
		para.status.Set_Current_Position(position);
		return(true);
	}
	else{
		para.status.Set_Show_Status(true);
//		ShowMessage(2,"no correct W", Jre_w, Jre_cut, re_mass);
		std::string position="no correct W  remass: "+Float_to_String(re_mass);
		para.status.Set_Current_Position(position);
		return(false);
	}
}

bool AN_RC_w_lnu(CDraw &para, float wmass, float mw_d, float mw_u, TLorentzVector &Vlep, TLorentzVector &Vmet, TLorentzVector &re_w, TLorentzVector &re_nu){
	para.status.Set_Current_Position("begin AN_RC_w_lnu");
	int lnum,mnum,key;


	TLorentzVector pw_choosed, pnu_choosed;



	float re_mass=0;
	//float xi_min=10000;
	bool Jre_w=false;

/// ShowMessage(2,"***************");
/// ShowMessage(2,"begin re W",wmass, mw_d, mw_u);
///	ShowMessage(2,"momentum",pneutrino, Vlep);
    int count=0;
	for(float wmass_tmp=0;wmass_tmp<5*wmass;wmass_tmp=wmass_tmp+1){
		TLorentzVector pw,  pnu;
		Jre_w=AN_RC_lnu(para, wmass_tmp, Vlep, Vmet, pw, pnu);
////    ShowMessage(2,"pnu",pnu);
////    ShowMessage(2,"pw",pw);
  		if(!Jre_w){
			continue;
  		}
		else{
	    	re_mass=pw.M();
	    	pw_choosed=pw;
	    	pnu_choosed=pnu;
			break;
		}
	////float xi=std::abs(pw.M()-wmass);
	////if(xi<xi_min){
	////	xi_min=xi;
	////	re_mass=pw.M();
	////	pw_choosed=pw;
	////	pnu_choosed=pnu;
	////}
	////else{
	////	count++;
	////}
	////ShowMessage(2,"wmass_tmp",wmass_tmp, xi, xi_min);
	////if(count>10){
	////	break;
	////}
	////std::string position="xi: "+Float_to_String(xi)+" xi_min: "+ Float_to_String(xi_min) + " re_mass: " + Float_to_String(re_mass)+ " pw.M: " + Float_to_String(pw.M());
	////para.status.Set_Current_Position(position);
  	}

	bool Jre_cut=false;

	if(re_mass>mw_d&&re_mass<mw_u){
		Jre_cut=true;
	}
	if(Jre_w&&Jre_cut){
		re_w .SetPxPyPzE(pw_choosed.Px(),pw_choosed.Py(),pw_choosed.Pz(),pw_choosed.E());
		re_nu.SetPxPyPzE(pnu_choosed.Px(),pnu_choosed.Py(),pnu_choosed.Pz(),pnu_choosed.E());
		if(re_w.M()<10){
			para.status.Set_Show_Status(true);
		}
		std::string position="has correct W  remass: "+Float_to_String(re_mass);
		para.status.Set_Current_Position(position);
		return(true);
	}
	else{
		para.status.Set_Show_Status(true);
//		ShowMessage(2,"no correct W", Jre_w, Jre_cut, re_mass);
		std::string position="no correct W  remass: "+Float_to_String(re_mass);
		para.status.Set_Current_Position(position);
		return(false);
	}
}

bool AN_RC_lnu(CDraw &para, float re_mass, TLorentzVector &plepton, TLorentzVector &pneutrino, TLorentzVector &pw, TLorentzVector &pnu){
	std::string position="begin AN_RC_lnu:"+Float_to_String(re_mass);
	para.status.Set_Current_Position(position);
	TLorentzVector pneu_plus = pneutrino;
	TLorentzVector pneu_minus= pneutrino;

	float neu_t = pneutrino.Pt();
	float lep_t = plepton.Pt();
	float lep_l = plepton.Pz();
	float lep_e = plepton.E();

	float var_alpha = std::pow(re_mass,2)+2*(plepton[0]*pneutrino[0]+plepton[1]*pneutrino[1]);
	float var_a = std::pow(lep_l,2)-std::pow(lep_e,2);
	float var_b = var_alpha*lep_l;
	float var_c = (std::pow(var_alpha,2) - 4*std::pow(lep_e,2)*std::pow(neu_t,2))/4;
	float var_d = std::pow(var_b,2) - 4*var_a*var_c;


	if(var_d < 0.0){
		return(false);
	}

	float neu_l_plus  = (-var_b+std::sqrt(var_d))/(2*var_a);
	float neu_l_minus = (-var_b-std::sqrt(var_d))/(2*var_a);

	pneu_plus (2) = neu_l_plus;
	pneu_plus (3) = std::sqrt(std::pow(neu_t,2)+std::pow(neu_l_plus,2));

	pneu_minus(2) = neu_l_minus;
	pneu_minus(3) = std::sqrt(std::pow(neu_t,2)+std::pow(neu_l_minus,2));

	TLorentzVector pw_plus  = plepton+pneu_plus;
	TLorentzVector pw_minus = plepton+pneu_minus;


	if(std::abs(pw_plus.M()-re_mass) < std::abs(pw_minus.M()-re_mass) ){
		pw =pw_plus;
		pnu=pneu_plus;
	} 
	else{ 
		pw =pw_minus;
		pnu=pneu_minus;
	}

//  ShowMessage(2,"remass",re_mass);
//  ShowMessage(2,"plepton",plepton);
//  ShowMessage(2,"pneutrino",pneutrino);
//  ShowMessage(2,"pneu_plus",pneu_plus);
//  ShowMessage(2,"pneu_minus",pneu_minus);
//  ShowMessage(2,"pw_plus",pw_plus);
//  ShowMessage(2,"pw_minus",pw_minus);
//  ShowMessage(2,"final pw ", pw);
//  ShowMessage(2,"final pnu ", pnu);
//  
	return(true);

}
